#Librerias necesarias 
library(astsa)
#install.packages(astsa)
#install.packages("ncdf4")
library(ncdf4)
#install.packages("data.table")
library(data.table)
#install.packages("forecast")
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## 
## Attaching package: 'forecast'
## The following object is masked from 'package:astsa':
## 
##     gas
#install.packages("urca")
library(urca)
#install.packages("rugarch")
library(rugarch)
## Loading required package: parallel
## 
## Attaching package: 'rugarch'
## The following object is masked from 'package:stats':
## 
##     sigma
#install.packages("plotly")
library(plotly)
## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout

INTRODUCCIÓN

Modelar la temperatura del aire puede ofrecer una comprensión más profunda de los patrones climáticos, permitiéndo hacer predicciones más precisas y contribuir al entendimiento de los cambios climáticos a largo plazo y sus implicaciones en el ámbito de la meteorología, la climatología y otras disciplinas relacionadas.

El meridiano cero, al estar en Greenwich, longitud cero (polo norte) es un lugar de referencia para la hora estándar mundial (GMT/UTC). La recopilación de datos meteorológicos en este lugar proporciona información clave para la creación y mantenimiento de registros climatológicos globales, comprender los patrones climáticos globales y las variaciones estacionales. Estos datos son fundamentales para la investigación sobre el cambio climático y sus efectos.

La temperatura del aire en el meridiano cero es utilizada en modelos meteorológicos globales para prever patrones climáticos y eventos meteorológicos extremos en todo el mundo. Utilizar el meridiano cero como referencia para mediciones atmosféricas proporciona estándares internacionales consistentes para la comparación de datos meteorológicos en todo el mundo e incluso puede ser crucial para la aviación y la navegación, ya que afecta a la densidad del aire, lo que a su vez influye en la altitud de vuelo y en el rendimiento de las aeronaves.

La medición de la temperatura del aire a una atmósfera de presión en el meridiano cero también puede contribuir a estudios atmosféricos más amplios, como la comprensión de la estratificación térmica y la circulación atmosférica.

La National Oceanic and Atmospheric Administration (NOAA) posee satélites para capturar mediciones de variables utilizando la técnica de asimilación de datos que combina observaciones del mundo real con modelos numéricos para mejorar la precisión y la confiabilidad de las predicciones o estimaciones.

En ese sentido se obtuvo una serie de tiempo con datos mensuales de temperatura promedio del aire en escala Kelvin a 1000 milibares de presión que equivalen a aproximadamente una atmósfera de presión (sobre el nivel del mar) durante los últimos 44 años (Enero de 1979 - Diciembre de 2022).

FUENTE: (https://psl.noaa.gov/data/gridded/data.ncep.reanalysis2.html):

#cargue de datos

ruta_archivo <-"/Users/macbookair/Documents/DOCTORADO/SERIES DE TIEMPO/PARCIAL 2/air.mon.mean.nc"
Vwind <- nc_open(ruta_archivo)

#En caso de error con los datos descargue los datos desde el sgte enlace y reemplace la ruta de archivo: https://drive.google.com/file/d/1PGXFUVWz9pC9_rCmmmfuJgjxe0N9j65L/view?usp=sharing

attributes(Vwind$var)
## $names
## [1] "time_bnds" "air"
attributes(Vwind$dim)
## $names
## [1] "lon"   "lat"   "level" "nbnds" "time"

Se extraen los datos de temperatura del aire para lon:1,lat:1 y level:1: Temperatura promedio del aire en escala Kelvin a 1000 milibares de presión (que equivalen a aproximadamente una atmósfera de presión o sobre el nivel del mar), en la ubicación geografica de latitud cero sobre el meridiano de Greenwich.

mi_variable <-ncvar_get(Vwind, "air")
nc_close(Vwind)
Vwind.ts <- ts(mi_variable[1,1,1,c(1:537)], frequency = 12)
str(Vwind.ts)
##  Time-Series [1:537] from 1 to 45.7: 245 239 248 251 263 ...
str(var)
## function (x, y = NULL, na.rm = FALSE, use)

Aplicando la metodología BOX-JENKINGS

IDENTIFICACIÓN: EDA

plot(decompose(Vwind.ts))

ggseasonplot(Vwind.ts)

ggseasonplot(Vwind.ts,polar=TRUE)

ggtsdisplay(Vwind.ts)

Modelo1<-sarima(Vwind.ts, p=1, d=1, q=1, P=1, D=1, Q=1, S=12)
## initial  value 1.363083 
## iter   2 value 1.052293
## iter   3 value 0.992212
## iter   4 value 0.913991
## iter   5 value 0.903168
## iter   6 value 0.895102
## iter   7 value 0.893648
## iter   8 value 0.892104
## iter   9 value 0.892088
## iter  10 value 0.892035
## iter  11 value 0.892032
## iter  12 value 0.892032
## iter  13 value 0.892031
## iter  14 value 0.892031
## iter  15 value 0.892031
## iter  15 value 0.892031
## iter  15 value 0.892031
## final  value 0.892031 
## converged
## initial  value 0.876949 
## iter   2 value 0.869482
## iter   3 value 0.859587
## iter   4 value 0.858531
## iter   5 value 0.857976
## iter   6 value 0.857340
## iter   7 value 0.857290
## iter   8 value 0.857278
## iter   9 value 0.857277
## iter  10 value 0.857277
## iter  10 value 0.857277
## final  value 0.857277 
## converged

Modelo1$ttable
##      Estimate     SE  t.value p.value
## ar1    0.2040 0.0494   4.1340  0.0000
## ma1   -0.9687 0.0231 -41.9310  0.0000
## sar1   0.0620 0.0513   1.2077  0.2277
## sma1  -0.9140 0.0274 -33.4061  0.0000
residuos<-Modelo1$fit$residuals
acf2(residuos)

##       [,1] [,2]  [,3] [,4]  [,5] [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13]
## ACF  -0.01 0.06 -0.01 0.05 -0.03 0.03 -0.01 -0.02 -0.03  0.09  0.01  0.00  0.04
## PACF -0.01 0.06 -0.01 0.05 -0.03 0.03 -0.01 -0.03 -0.03  0.09  0.02 -0.01  0.04
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF   0.00  0.03  0.02  0.00 -0.03  0.02  0.03 -0.02     0  0.04 -0.01  0.01
## PACF -0.01  0.03  0.01 -0.01 -0.03  0.02  0.02 -0.02     0  0.03 -0.01  0.00
##      [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37]
## ACF  -0.01 -0.04 -0.08 -0.01 -0.02  0.00 -0.04  0.05  0.07 -0.05  0.07 -0.01
## PACF -0.02 -0.05 -0.07 -0.01 -0.02  0.01 -0.04  0.05  0.09 -0.07  0.06  0.00
##      [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## ACF  -0.01 -0.02  0.02 -0.02  0.02  0.02 -0.06 -0.02  0.00  0.07 -0.02
## PACF -0.01 -0.01  0.02  0.00  0.03  0.03 -0.08  0.00 -0.01  0.08 -0.01
Box.test(residuos, lag=8, type='Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  residuos
## X-squared = 4.9717, df = 8, p-value = 0.7606
Box.test(residuos, lag=16, type='Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  residuos
## X-squared = 11.285, df = 16, p-value = 0.7915
Box.test(residuos, lag=24, type='Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  residuos
## X-squared = 13.441, df = 24, p-value = 0.9583

Debo verificar si hay raices unitarias en las estacionalidades y volatilidad. De la exploración de los datos, en las diferentes gráficas se puede identificar claramente un componente estacional en los datos originales por lo que una vez verificada la estacionariedad, el modelo propuesto debe poder representar la estacionalidad.

ESTIMACIÓN: Detección de raíces unitarias y verificar estabilidad

ESTACIONARIEDAD

Se verifica estacionariedad para serie original siguiendo el algoritmo de la prueba de Dickey Fuller aumentada:

ADF<-ur.df(Vwind.ts,type = "trend",lag=12)
summary(ADF)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.3053 -1.2500  0.0843  1.4492  9.5314 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  93.202007  29.494369   3.160  0.00167 ** 
## z.lag.1      -0.363249   0.115071  -3.157  0.00169 ** 
## tt            0.003510   0.001379   2.546  0.01120 *  
## z.diff.lag1  -0.300758   0.116838  -2.574  0.01033 *  
## z.diff.lag2  -0.361666   0.111816  -3.234  0.00130 ** 
## z.diff.lag3  -0.470242   0.103546  -4.541 6.98e-06 ***
## z.diff.lag4  -0.455917   0.093954  -4.853 1.62e-06 ***
## z.diff.lag5  -0.528129   0.085185  -6.200 1.17e-09 ***
## z.diff.lag6  -0.527841   0.077591  -6.803 2.89e-11 ***
## z.diff.lag7  -0.557240   0.070112  -7.948 1.23e-14 ***
## z.diff.lag8  -0.619246   0.062223  -9.952  < 2e-16 ***
## z.diff.lag9  -0.692397   0.056998 -12.148  < 2e-16 ***
## z.diff.lag10 -0.622403   0.052706 -11.809  < 2e-16 ***
## z.diff.lag11 -0.422243   0.049733  -8.490 2.27e-16 ***
## z.diff.lag12 -0.080394   0.043949  -1.829  0.06795 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.501 on 509 degrees of freedom
## Multiple R-squared:  0.8399, Adjusted R-squared:  0.8355 
## F-statistic: 190.8 on 14 and 509 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -3.1567 3.4028 5.0195 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.96 -3.41 -3.12
## phi2  6.09  4.68  4.03
## phi3  8.27  6.25  5.34

Para evaluar la presencia de raíces unitarias en una serie temporal, lo que indica la no estacionariedad se seguirán los sigueintes pasos: Paso 1. Iniciar con el modelo menos restrictivo con un modelo que incluye tendencia,, usar \({\tau_3}\) para probar \({\gamma=0}\). Dado que la prueba tiene poca potencia si se rechaza entonces no se continua.

ADF1<-ur.df(Vwind.ts,type = "trend",lag=24, selectlags= "BIC")
summary(ADF1)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.3048 -1.3445  0.1275  1.3024  9.4597 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  106.202015  29.879545   3.554 0.000415 ***
## z.lag.1       -0.414247   0.116625  -3.552 0.000419 ***
## tt             0.004223   0.001440   2.932 0.003519 ** 
## z.diff.lag1   -0.221842   0.114052  -1.945 0.052325 .  
## z.diff.lag2   -0.271234   0.105311  -2.576 0.010296 *  
## z.diff.lag3   -0.372543   0.094459  -3.944 9.16e-05 ***
## z.diff.lag4   -0.369191   0.085311  -4.328 1.82e-05 ***
## z.diff.lag5   -0.448534   0.076825  -5.838 9.51e-09 ***
## z.diff.lag6   -0.449656   0.068732  -6.542 1.51e-10 ***
## z.diff.lag7   -0.481309   0.059761  -8.054 5.96e-15 ***
## z.diff.lag8   -0.547546   0.052546 -10.420  < 2e-16 ***
## z.diff.lag9   -0.630295   0.045861 -13.743  < 2e-16 ***
## z.diff.lag10  -0.567666   0.044048 -12.887  < 2e-16 ***
## z.diff.lag11  -0.370945   0.041787  -8.877  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.523 on 498 degrees of freedom
## Multiple R-squared:  0.8374, Adjusted R-squared:  0.8331 
## F-statistic: 197.3 on 13 and 498 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -3.552 4.2529 6.3239 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.96 -3.41 -3.12
## phi2  6.09  4.68  4.03
## phi3  8.27  6.25  5.34
ADF2<-ur.df(Vwind.ts,type = "trend",lag=24, selectlags= "AIC")
summary(ADF2)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.8099 -1.2605  0.0333  1.4008  8.9909 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  95.865549  34.377315   2.789 0.005501 ** 
## z.lag.1      -0.373916   0.134193  -2.786 0.005538 ** 
## tt            0.003898   0.001582   2.465 0.014056 *  
## z.diff.lag1  -0.323768   0.135091  -2.397 0.016922 *  
## z.diff.lag2  -0.330053   0.132635  -2.488 0.013164 *  
## z.diff.lag3  -0.394377   0.130510  -3.022 0.002645 ** 
## z.diff.lag4  -0.361216   0.129402  -2.791 0.005454 ** 
## z.diff.lag5  -0.407340   0.128011  -3.182 0.001556 ** 
## z.diff.lag6  -0.386651   0.126576  -3.055 0.002377 ** 
## z.diff.lag7  -0.418297   0.125470  -3.334 0.000922 ***
## z.diff.lag8  -0.460856   0.124773  -3.694 0.000246 ***
## z.diff.lag9  -0.518100   0.124557  -4.160 3.77e-05 ***
## z.diff.lag10 -0.454651   0.124858  -3.641 0.000300 ***
## z.diff.lag11 -0.349479   0.125576  -2.783 0.005595 ** 
## z.diff.lag12 -0.110921   0.124937  -0.888 0.375080    
## z.diff.lag13 -0.031793   0.120786  -0.263 0.792492    
## z.diff.lag14 -0.136572   0.114844  -1.189 0.234943    
## z.diff.lag15 -0.142753   0.107624  -1.326 0.185330    
## z.diff.lag16 -0.166168   0.100754  -1.649 0.099742 .  
## z.diff.lag17 -0.184907   0.093931  -1.969 0.049573 *  
## z.diff.lag18 -0.232291   0.087293  -2.661 0.008048 ** 
## z.diff.lag19 -0.212242   0.079592  -2.667 0.007917 ** 
## z.diff.lag20 -0.196943   0.071869  -2.740 0.006364 ** 
## z.diff.lag21 -0.254496   0.062229  -4.090 5.06e-05 ***
## z.diff.lag22 -0.292448   0.054030  -5.413 9.77e-08 ***
## z.diff.lag23 -0.161028   0.044749  -3.598 0.000353 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.438 on 486 degrees of freedom
## Multiple R-squared:  0.8519, Adjusted R-squared:  0.8443 
## F-statistic: 111.8 on 25 and 486 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -2.7864 2.7515 3.8821 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.96 -3.41 -3.12
## phi2  6.09  4.68  4.03
## phi3  8.27  6.25  5.34

Teniendo en cuenta queValor del estadistico de prueba -2.7864 mayor que tau3 critico para cualquier nivel de significancia, por tanto no se rechaza Ho

Se selecciona el ajuste de ADF2 escogido por AIC con 23 rezagos o grados de libertad, cuyos residuos se comportan más como ruido blanco comparado con el ajuste ADF1 escogido por BIC con 11 rezagos. El summary de este ajuste arroja un valor del estadístico de -2.7864 el cual comparado con tau3 critico para cualquier nivel de significancia no permite rechazar la hipótesis nula. Este no es el resultado final ya que el proposito de esta parte es determinar el numero de lags. El pvalor dado en la tabla t para los coeficientes en particular el de gamma no es confiable porque la distribucion de los procesos con raiz unitaria no se comporta como t.

p-valor dado en la tabla t para los coeficientes, en particular el de gamma, no es confiable debido a que la distribución de los procesos con raíz unitaria no se comporta como una distribución t. La prueba de Dickey-Fuller es especialmente sensible a la cantidad de rezagos utilizados, y seleccionar el número óptimo de rezagos es una parte importante del proceso.

Se debe probar el modelo nuevamente según el paso 2 de la imagen para determinar si la serie no es estacionaria.

residuos<-ADF1@res
acf2(residuos)

##       [,1] [,2] [,3] [,4] [,5] [,6] [,7]  [,8] [,9] [,10] [,11] [,12] [,13]
## ACF  -0.03 0.02 0.04 0.05 0.03 0.03    0  0.00 0.01  0.01 -0.08 -0.18  0.00
## PACF -0.03 0.02 0.04 0.05 0.04 0.03    0 -0.01 0.00  0.01 -0.08 -0.19 -0.02
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF  -0.06  0.05  0.06  0.03 -0.02  0.03  0.06 -0.04 -0.09  0.02  0.00  0.00
## PACF -0.05  0.07  0.10  0.07  0.00  0.03  0.05 -0.05 -0.12 -0.04 -0.05 -0.02
##      [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## ACF  -0.06 -0.06 -0.06 -0.01 -0.02  0.01 -0.01  0.05
## PACF -0.07 -0.01 -0.01  0.03  0.01  0.05  0.02  0.03
residuos<-ADF2@res
acf2(residuos)

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]  [,9] [,10] [,11] [,12] [,13]
## ACF     0 0.01 0.04 0.03 0.01 0.01 0.01    0 -0.01 -0.02 -0.02 -0.07 -0.06
## PACF    0 0.01 0.04 0.03 0.01 0.00 0.01    0 -0.01 -0.02 -0.03 -0.07 -0.06
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF   0.02  0.07  0.06  0.02  0.03  0.02  0.03  0.01  0.01 -0.08 -0.14 -0.04
## PACF  0.03  0.08  0.07  0.02  0.02  0.01  0.02  0.00  0.00 -0.10 -0.16 -0.05
##      [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## ACF  -0.06 -0.04 -0.05 -0.01 -0.01  0.04 -0.01  0.06
## PACF -0.05 -0.01 -0.02  0.01  0.00  0.05  0.01  0.07

Se visualiza la autocorrelación de los residuos y evaluar su comportamiento.

# Hacemos la prueba de Dickey Fuller con el modelo restringido.

ADF_R<-ur.df(Vwind.ts,type = "drift",lag=23, selectlags= "Fixed")
summary(ADF_R)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.747 -1.284  0.025  1.325  8.889 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  20.96201   16.08693   1.303 0.193173    
## z.lag.1      -0.08061    0.06209  -1.298 0.194784    
## z.diff.lag1  -0.60300    0.07376  -8.175 2.57e-15 ***
## z.diff.lag2  -0.59752    0.07653  -7.807 3.60e-14 ***
## z.diff.lag3  -0.64994    0.07937  -8.188 2.33e-15 ***
## z.diff.lag4  -0.60562    0.08337  -7.264 1.49e-12 ***
## z.diff.lag5  -0.64080    0.08644  -7.413 5.49e-13 ***
## z.diff.lag6  -0.60744    0.08977  -6.767 3.79e-11 ***
## z.diff.lag7  -0.62689    0.09297  -6.743 4.40e-11 ***
## z.diff.lag8  -0.65671    0.09651  -6.805 2.97e-11 ***
## z.diff.lag9  -0.70132    0.10034  -6.989 9.12e-12 ***
## z.diff.lag10 -0.62376    0.10465  -5.961 4.82e-09 ***
## z.diff.lag11 -0.50635    0.10859  -4.663 4.03e-06 ***
## z.diff.lag12 -0.25580    0.11032  -2.319 0.020822 *  
## z.diff.lag13 -0.16559    0.10803  -1.533 0.125983    
## z.diff.lag14 -0.25922    0.10368  -2.500 0.012742 *  
## z.diff.lag15 -0.25192    0.09836  -2.561 0.010733 *  
## z.diff.lag16 -0.26199    0.09316  -2.812 0.005118 ** 
## z.diff.lag17 -0.26794    0.08773  -3.054 0.002380 ** 
## z.diff.lag18 -0.30302    0.08249  -3.673 0.000266 ***
## z.diff.lag19 -0.27037    0.07605  -3.555 0.000414 ***
## z.diff.lag20 -0.24380    0.06927  -3.520 0.000472 ***
## z.diff.lag21 -0.29010    0.06063  -4.784 2.28e-06 ***
## z.diff.lag22 -0.31674    0.05303  -5.972 4.51e-09 ***
## z.diff.lag23 -0.17469    0.04459  -3.918 0.000102 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.448 on 488 degrees of freedom
## Multiple R-squared:   0.85,  Adjusted R-squared:  0.8426 
## F-statistic: 115.2 on 24 and 488 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -1.2983 1.0894 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1  6.43  4.59  3.78

Paso 2. Verificar si debe incluirse la tendencia, usando \({\phi_3}\). Probemos entonces que \({a_2 = \gamma = 0}\). Recordemos las hipótesis involucradas y el cálculo de \({\phi_i}\)

SSR_U<-ADF2@testreg$sigma^2*ADF2@testreg$df[2]
SSR_R<-ADF_R@testreg$sigma^2*ADF_R@testreg$df[2]
r<-2 #Numero de restricciones
df_u<-ADF@testreg$df[2]
fi3=((SSR_R-SSR_U)/r)/(SSR_U/df_u)
fi3
## [1] 3.182458

Se calcula SSR para los modelos restringido y no restringido, define el número de restricciones y los grados de libertad, y luego calcula el estadístico \({\phi_3}\). con el fin de evaluar si la inclusión de la tendencia mejora significativamente el modelo, comparando el modelo no restringido (con tendencia) y el modelo restringido (sin tendencia).

Se espera que \({\phi_3}\) sea mayor a los puntos críticos del modelo no restringido para rechazar la hipótesis nula, se obtiene que es \({\phi_3}\) 3.038654 Al comparar con los puntos críticos del modelo ADF2 (NO RESTRINGIDO) se tiene que \({\phi_3}\) es menor por lo que por tanto no se rechaza H0, Quiere decir que entonces que no hay tendencia, puede haber raíz unitaria y debemos ir al siguiente paso del algoritmo.

# Miramos la ACF de los residuos para estar seguro que hay ruido blanco
residuos_ADF_R<-ADF_R@res
acf2(residuos_ADF_R)

##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]  [,9] [,10] [,11] [,12] [,13]
## ACF  -0.01 0.01 0.03 0.03 0.01 0.01 0.01 0.01 -0.01 -0.02 -0.02 -0.06 -0.06
## PACF -0.01 0.01 0.03 0.03 0.01 0.00 0.01 0.01 -0.01 -0.02 -0.02 -0.07 -0.06
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF   0.02  0.07  0.05  0.01  0.02  0.01  0.03     0  0.01 -0.08 -0.15 -0.04
## PACF  0.02  0.07  0.06  0.02  0.02  0.01  0.02     0  0.00 -0.10 -0.17 -0.05
##      [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## ACF  -0.06 -0.04 -0.05 -0.01 -0.01  0.04 -0.01  0.06
## PACF -0.05 -0.01 -0.02  0.01  0.01  0.06  0.01  0.07

Al revisar la ACF de los residuos vemos que hay una espiga que sobresale un poco y puede ser producto del azar. Por lo tanto, asumiremos que la ACF de los residuos es ruido blanco.

Paso 3. Estimaremos el modelo sin tendencia, es decir, \({\Delta Y_t = a_0 + \gamma * Y_{t-1}+c_1* \Delta Y{t-1} + ... + cp * \Delta Y_{t-p} + \epsilon_t}\). Para eso usamos \({\tau2}\) para verificar si hay raiz unitaria.

summary(ADF_R)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.747 -1.284  0.025  1.325  8.889 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  20.96201   16.08693   1.303 0.193173    
## z.lag.1      -0.08061    0.06209  -1.298 0.194784    
## z.diff.lag1  -0.60300    0.07376  -8.175 2.57e-15 ***
## z.diff.lag2  -0.59752    0.07653  -7.807 3.60e-14 ***
## z.diff.lag3  -0.64994    0.07937  -8.188 2.33e-15 ***
## z.diff.lag4  -0.60562    0.08337  -7.264 1.49e-12 ***
## z.diff.lag5  -0.64080    0.08644  -7.413 5.49e-13 ***
## z.diff.lag6  -0.60744    0.08977  -6.767 3.79e-11 ***
## z.diff.lag7  -0.62689    0.09297  -6.743 4.40e-11 ***
## z.diff.lag8  -0.65671    0.09651  -6.805 2.97e-11 ***
## z.diff.lag9  -0.70132    0.10034  -6.989 9.12e-12 ***
## z.diff.lag10 -0.62376    0.10465  -5.961 4.82e-09 ***
## z.diff.lag11 -0.50635    0.10859  -4.663 4.03e-06 ***
## z.diff.lag12 -0.25580    0.11032  -2.319 0.020822 *  
## z.diff.lag13 -0.16559    0.10803  -1.533 0.125983    
## z.diff.lag14 -0.25922    0.10368  -2.500 0.012742 *  
## z.diff.lag15 -0.25192    0.09836  -2.561 0.010733 *  
## z.diff.lag16 -0.26199    0.09316  -2.812 0.005118 ** 
## z.diff.lag17 -0.26794    0.08773  -3.054 0.002380 ** 
## z.diff.lag18 -0.30302    0.08249  -3.673 0.000266 ***
## z.diff.lag19 -0.27037    0.07605  -3.555 0.000414 ***
## z.diff.lag20 -0.24380    0.06927  -3.520 0.000472 ***
## z.diff.lag21 -0.29010    0.06063  -4.784 2.28e-06 ***
## z.diff.lag22 -0.31674    0.05303  -5.972 4.51e-09 ***
## z.diff.lag23 -0.17469    0.04459  -3.918 0.000102 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.448 on 488 degrees of freedom
## Multiple R-squared:   0.85,  Adjusted R-squared:  0.8426 
## F-statistic: 115.2 on 24 and 488 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -1.2983 1.0894 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1  6.43  4.59  3.78

Al comparar tau_2 se tiene que es -1,2983 con los valores críticos de tau2 y se tiene que patra para cualquier nivel de significancia es mayor por lo que no se rechaza la hipótesis nula. Lo que indica que aún es posible que exista una raíz unitaria. y debemos calcular phi1.

ADF_N<-ur.df(Vwind.ts,type = "none",lag=23, selectlags= "Fixed")
summary(ADF_N)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.8253 -1.2291  0.0453  1.3413  8.6562 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## z.lag.1       0.0002904  0.0004190   0.693 0.488653    
## z.diff.lag1  -0.6796656  0.0445179 -15.267  < 2e-16 ***
## z.diff.lag2  -0.6706184  0.0520994 -12.872  < 2e-16 ***
## z.diff.lag3  -0.7194790  0.0587964 -12.237  < 2e-16 ***
## z.diff.lag4  -0.6716356  0.0662581 -10.137  < 2e-16 ***
## z.diff.lag5  -0.7034662  0.0718751  -9.787  < 2e-16 ***
## z.diff.lag6  -0.6667707  0.0774220  -8.612  < 2e-16 ***
## z.diff.lag7  -0.6832754  0.0823424  -8.298 1.04e-15 ***
## z.diff.lag8  -0.7101975  0.0874058  -8.125 3.68e-15 ***
## z.diff.lag9  -0.7519159  0.0925904  -8.121 3.80e-15 ***
## z.diff.lag10 -0.6712112  0.0981773  -6.837 2.42e-11 ***
## z.diff.lag11 -0.5514277  0.1030060  -5.353 1.33e-07 ***
## z.diff.lag12 -0.2982144  0.1054810  -2.827 0.004888 ** 
## z.diff.lag13 -0.2053551  0.1037045  -1.980 0.048242 *  
## z.diff.lag14 -0.2961732  0.0997991  -2.968 0.003148 ** 
## z.diff.lag15 -0.2855764  0.0949771  -3.007 0.002776 ** 
## z.diff.lag16 -0.2923978  0.0902530  -3.240 0.001278 ** 
## z.diff.lag17 -0.2948816  0.0853173  -3.456 0.000595 ***
## z.diff.lag18 -0.3266302  0.0805314  -4.056 5.81e-05 ***
## z.diff.lag19 -0.2901430  0.0745697  -3.891 0.000114 ***
## z.diff.lag20 -0.2597082  0.0682293  -3.806 0.000159 ***
## z.diff.lag21 -0.3022332  0.0599564  -5.041 6.54e-07 ***
## z.diff.lag22 -0.3249980  0.0526919  -6.168 1.45e-09 ***
## z.diff.lag23 -0.1791509  0.0444872  -4.027 6.55e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.449 on 489 degrees of freedom
## Multiple R-squared:  0.8495, Adjusted R-squared:  0.8421 
## F-statistic:   115 on 24 and 489 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: 0.693 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62
SSR_U<-ADF_R@testreg$sigma^2*ADF_R@testreg$df[2]
SSR_R<-ADF_N@testreg$sigma^2*ADF_N@testreg$df[2]
r<-2 #Numero de restricciones
df_u<-ADF_R@testreg$df[2]
fi1=((SSR_R-SSR_U)/r)/(SSR_U/df_u)
fi1
## [1] 0.8489648
test_001<-fi1>ADF_R@cval[2,1]
test_005<-fi1>ADF_R@cval[2,2]
test_010<-fi1>ADF_R@cval[2,3]
# Miramos la ACF de los residuos para estar seguro que hay ruido blanco
residuos_ADF_N<-ADF_N@res
acf2(residuos_ADF_N)

##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]  [,9] [,10] [,11] [,12] [,13]
## ACF  -0.01 0.01 0.03 0.03 0.01    0 0.01 0.01 -0.01 -0.02 -0.02 -0.07 -0.06
## PACF -0.01 0.01 0.03 0.03 0.01    0 0.01 0.00 -0.01 -0.02 -0.02 -0.07 -0.06
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF   0.01  0.06  0.04  0.01  0.02  0.01  0.02  0.00  0.00 -0.09 -0.16 -0.04
## PACF  0.01  0.07  0.06  0.01  0.01  0.00  0.02 -0.01 -0.01 -0.10 -0.17 -0.06
##      [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## ACF  -0.06 -0.04 -0.06 -0.02 -0.01  0.04 -0.01  0.06
## PACF -0.06 -0.01 -0.03  0.00  0.00  0.05  0.00  0.07

Teniendo en cuenta que phi_1 es menor para cualquier nivel de confianza que los puntos críticos del modelo no restringido para este punto, y los residuos se comportan como ruido blanco, entonces se acepta la hipótesis nula y se concluye que la serie de tiempo tiene tendencia estocástica por lo que hay raiz unitaria, sin intersecto ni tendencia.

DVwind<-(diff(Vwind.ts))
plot(DVwind)

Viendo el resultado que arroja la gráfica de la serie con una diferenciación, vemos un comportamiento que fluctua alrededor de 0.Ahora bien, para tener mayor certeza que la serie una vez diferenciada es estacionaria se aplica la función ndiffs que permite determinar el número de diferenciaciones que requiere la serie para volverse estacionaria, el resultado es 1:

ndiffs(Vwind.ts,alpha = 0.05,test = c("kpss", "adf", "pp"),type = c("level", "trend"),max.d = 2)
## [1] 1

Se revisan los residuos:

acf2(DVwind)

##      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12]
## ACF  0.55  0.28 -0.09 -0.35 -0.57 -0.61 -0.56 -0.37 -0.09  0.30  0.61  0.78
## PACF 0.55 -0.02 -0.34 -0.27 -0.32 -0.29 -0.35 -0.39 -0.44 -0.33 -0.17  0.16
##      [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF   0.62  0.27 -0.07 -0.36 -0.54 -0.62 -0.54 -0.35 -0.08  0.27  0.62  0.76
## PACF  0.20  0.00 -0.03 -0.04 -0.01 -0.05  0.01  0.03 -0.06 -0.15 -0.03  0.09
##      [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
## ACF   0.62  0.28 -0.07 -0.36 -0.53 -0.60 -0.53 -0.37 -0.06  0.28  0.58  0.76
## PACF  0.10  0.01 -0.04 -0.06  0.01  0.01  0.03 -0.03  0.00 -0.01 -0.10  0.09
##      [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## ACF   0.60  0.27 -0.08 -0.34 -0.53 -0.58 -0.52 -0.36 -0.07  0.26  0.61  0.71
## PACF  0.05 -0.04 -0.06 -0.01 -0.02  0.01  0.06  0.01  0.02 -0.09  0.06  0.00

Para entender si hay patrones de correlación en los datos diferenciados se realiza el gráfico de autocorrelación de la serie diferenciada. En este se puede observar una autocorrelación Positiva en el Primer Rezago sugiere que hay una correlación positiva entre los valores actuales y los valores inmediatamente anteriores. En otras palabras, la temperatura del aire en un momento dado está correlacionada positivamente con la temperatura del aire en el momento anterior. Esto seguido de “ocho espigas hacia abajo” podrían indicar autocorrelaciones negativas en los rezagos siguientes. Esto podría interpretarse como una especie de patrón cíclico de cambios en la temperatura del aire. Ahora bien, los rezagos están por fuera de las bandas de confianza (las “bandas” generadas por la función acf()), podría indicar que algunas de estas autocorrelaciones son estadísticamente significativas.

Un patrón claro en las autocorrelaciones podría indica estacionalidad o repetición de patrones en la serie temporal de la temperatura del aire, lo que puede ser resultado de características específicas los datos de temperatura del aire, como factores climáticos, estacionales, o patrones diarios pueden influir en estos resultados. Se ratifica el comportamiento estacional de la serie aún despúes de una diferenciación y al no hallar ruido blanco implica que la diferenciación no hace que la serie sea estacionaria, en ese caso se se requiere probar una diferenciación estacional para cada 12 periodos:

plot(diff((Vwind.ts),lag=12))

ggseasonplot(diff(diff(Vwind.ts),lag=12))

ggseasonplot(diff(diff(Vwind.ts),lag=12), polar=TRUE)

ggsubseriesplot(diff(diff(Vwind.ts),lag=12))

acf2(diff(diff(Vwind.ts), lag=12))

##       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12]
## ACF  -0.43 -0.04 -0.06  0.07 -0.06  0.05  0.00 -0.05 -0.01  0.10  0.15 -0.46
## PACF -0.43 -0.28 -0.26 -0.14 -0.17 -0.10 -0.06 -0.12 -0.13  0.01  0.30 -0.29
##      [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF   0.22 -0.03  0.06  0.00  0.02 -0.07  0.02  0.06 -0.03 -0.08  0.11 -0.07
## PACF -0.14 -0.17 -0.14 -0.04 -0.05 -0.06 -0.04  0.02 -0.03 -0.04  0.27 -0.16
##      [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
## ACF   0.02  0.03 -0.01 -0.07  0.05  0.01 -0.01 -0.05  0.06  0.08 -0.16  0.14
## PACF -0.06 -0.06 -0.04 -0.04 -0.03 -0.05 -0.08 -0.06 -0.07  0.11  0.12 -0.03
##      [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## ACF  -0.08  0.00  0.00  0.05 -0.08  0.06  0.04 -0.04 -0.04 -0.03  0.11 -0.11
## PACF -0.04 -0.09 -0.04 -0.05 -0.10 -0.06  0.04 -0.02 -0.04  0.02  0.14 -0.02

Se hace la diferenciación estacional de 12 periodos al tratarse de datos mensuales

d12datos<-diff(Vwind.ts, lag=12)
acf2(d12datos)

##      [,1] [,2]  [,3] [,4]  [,5] [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13]
## ACF  0.18 0.05  0.01 0.05 -0.01 0.03 -0.01 -0.05 -0.01  0.04 -0.07 -0.43 -0.03
## PACF 0.18 0.02 -0.01 0.05 -0.03 0.03 -0.02 -0.05  0.02  0.04 -0.08 -0.42  0.14
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF  -0.01  0.08  0.06  0.04 -0.01  0.04  0.05 -0.03 -0.07  0.01 -0.08 -0.04
## PACF  0.02  0.09  0.08 -0.02  0.01  0.02  0.01 -0.06 -0.01  0.00 -0.33  0.08
##      [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37]
## ACF  -0.03 -0.08 -0.11 -0.02 -0.01 -0.02 -0.01  0.08  0.07 -0.07  0.05 -0.05
## PACF -0.01  0.01 -0.01  0.00 -0.01  0.01  0.05  0.03  0.05 -0.13 -0.16 -0.02
##      [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## ACF  -0.02  0.02  0.05  0.00  0.07  0.05 -0.04 -0.05 -0.01  0.08  0.00
## PACF -0.02  0.05  0.01  0.03  0.09  0.07 -0.04  0.02  0.05 -0.01 -0.13
muestra<-window(d12datos, end=c(430,12))
## Warning in window.default(x, ...): 'end' value not changed

Nótese que desaparecen las espigas estacionales. Se repite el proceso de estimación del número de rezagos en diferencia

TryMod<-ur.df(muestra, type="trend", lags=12, selectlags = "AIC")
summary(TryMod)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.9712  -1.3831   0.0122   1.4974   9.9926 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   9.775e-02  2.464e-01   0.397 0.691798    
## z.lag.1      -1.077e+00  1.287e-01  -8.368 5.96e-16 ***
## tt           -5.024e-05  8.044e-04  -0.062 0.950225    
## z.diff.lag1   2.674e-01  1.169e-01   2.287 0.022609 *  
## z.diff.lag2   3.243e-01  1.129e-01   2.874 0.004230 ** 
## z.diff.lag3   2.995e-01  1.092e-01   2.743 0.006310 ** 
## z.diff.lag4   3.295e-01  1.042e-01   3.162 0.001660 ** 
## z.diff.lag5   2.995e-01  9.899e-02   3.026 0.002607 ** 
## z.diff.lag6   3.473e-01  9.350e-02   3.714 0.000227 ***
## z.diff.lag7   3.230e-01  8.808e-02   3.667 0.000272 ***
## z.diff.lag8   2.883e-01  8.096e-02   3.562 0.000404 ***
## z.diff.lag9   2.761e-01  7.361e-02   3.752 0.000196 ***
## z.diff.lag10  3.474e-01  6.467e-02   5.372 1.20e-07 ***
## z.diff.lag11  3.273e-01  5.628e-02   5.815 1.09e-08 ***
## z.diff.lag12 -1.273e-01  4.475e-02  -2.845 0.004625 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.687 on 497 degrees of freedom
## Multiple R-squared:  0.5424, Adjusted R-squared:  0.5295 
## F-statistic: 42.08 on 14 and 497 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -8.3684 23.3628 35.0426 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.96 -3.41 -3.12
## phi2  6.09  4.68  4.03
## phi3  8.27  6.25  5.34
resTryMod<-TryMod@res
acf2(resTryMod, max.lag = 30)

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]  [,9] [,10] [,11] [,12] [,13]
## ACF     0    0 0.03 0.03    0    0 0.02 0.01 -0.03 -0.01  0.04 -0.17 -0.02
## PACF    0    0 0.03 0.03    0    0 0.02 0.01 -0.03 -0.01  0.04 -0.17 -0.01
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF   0.01  0.06  0.04     0  0.02  0.01  0.03 -0.04  0.03  0.00 -0.31     0
## PACF  0.01  0.07  0.05     0  0.01  0.00  0.04 -0.05  0.02  0.01 -0.35     0
##      [,26] [,27] [,28] [,29] [,30]
## ACF  -0.02 -0.03 -0.09     0  0.00
## PACF -0.02  0.01 -0.05     0  0.01
Box.test(resTryMod, lag=8, type='Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  resTryMod
## X-squared = 1.3205, df = 8, p-value = 0.9953
Box.test(resTryMod, lag=16, type='Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  resTryMod
## X-squared = 20.569, df = 16, p-value = 0.1957
Box.test(resTryMod, lag=24, type='Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  resTryMod
## X-squared = 74.678, df = 24, p-value = 4.186e-07
#Revisión de volatilidad
vol<-100*diff(log(Vwind.ts))
vol<-na.omit(vol)
plot(vol)

Una vez calculado y graficado los rendimientos porcentuales de la temperatura del aire, para obtener información sobre la volatilidad de la temperatura del aire a lo largo del tiempo. No se observan patrones notables en el gráfico, como picos o caídas bruscas, por lo que no existen indicios que puedan interpretarse como momentos de mayor o menor volatilidad en la temperatura del aire.

# Ajusta un modelo GARCH(1,1)
model_spec <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(0, 1)),
                         mean.model = list(armaOrder = c(0, 0)))

fit <- ugarchfit(spec = model_spec, data = d12datos)

# Obtiene los residuos estandarizados
residuals <- residuals(fit, standardize = TRUE)

# Crea un gráfico de dispersión condicional
plot(x = d12datos, y = residuals, type = "l",
     xlab = "Tiempo", ylab = "Residuos Estandarizados",
     main = "Gráfico de Dispersión Condicional")

# Añade una línea de referencia en y=0 para destacar los cambios en la volatilidad
abline(h = 0, col = "red", lty = 2)

Adicionalmente para analizar la volatilidad condicional en la serie de tiempo se realiza una gráfica de dispersión condicional que muestra la evolución de los residuos estandarizados a lo largo del tiempo que refleja que no hay mayor volatilidad debido a que gráficamente no se visualizan cambios destacados con respecto a la línea de referencia en y=0.

MODELACIÓN Y AJUSTE

A sabienda que existe estacionalidad en los datos originales, vamos a tratar de construir el modelo a partir de la misma, esto nos supone que en el modelo SARIMA \(d=1\).

Inicialmente se usa la función arima para ajustar un modelo ARIMA estacional. El orden del componente no estacional seleccionado como punto de partida según el análisis de autocorrelación de la serie de tiempo es (2,1,0) (indicando un AR(2) y una diferenciación de orden 1). El componente estacional se especifica como (0,0,5) (indicando un componente estacional de orden 5).

M1<-arima(muestra,order=c(2,1,0),seasonal=list(order=c(0,0,5)),fixed = c(0,NA,0,NA,0,0,NA))
## Warning in arima(muestra, order = c(2, 1, 0), seasonal = list(order = c(0, :
## some AR parameters were fixed: setting transform.pars = FALSE
M1
## 
## Call:
## arima(x = muestra, order = c(2, 1, 0), seasonal = list(order = c(0, 0, 5)), 
##     fixed = c(0, NA, 0, NA, 0, 0, NA))
## 
## Coefficients:
##       ar1      ar2  sma1     sma2  sma3  sma4    sma5
##         0  -0.0474     0  -0.1100     0     0  0.0537
## s.e.    0   0.0438     0   0.0493     0     0  0.0486
## 
## sigma^2 estimated as 14.98:  log likelihood = -1452.97,  aic = 2913.94
bic=AIC(M1,k = log(length(d12datos)))
bic
## [1] 2930.994
acf2(M1$residuals)

##       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12]
## ACF  -0.44  0.00 -0.09  0.06 -0.05  0.05 -0.01 -0.04  0.01  0.08  0.17 -0.50
## PACF -0.44 -0.24 -0.25 -0.15 -0.17 -0.10 -0.07 -0.12 -0.11  0.02  0.34 -0.31
##      [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF   0.23 -0.04  0.08  0.00  0.02 -0.06  0.03  0.05 -0.03 -0.08  0.07  0.01
## PACF -0.18 -0.17 -0.14 -0.04 -0.06 -0.06 -0.03  0.02 -0.02 -0.05  0.23 -0.14
##      [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
## ACF   0.00  0.02 -0.02 -0.07  0.05  0.00 -0.01 -0.04  0.05  0.09 -0.14  0.09
## PACF -0.06 -0.05 -0.03 -0.04 -0.02 -0.05 -0.08 -0.04 -0.05  0.10  0.10 -0.02
##      [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## ACF  -0.07  0.00  0.01  0.06 -0.08  0.06  0.04 -0.04 -0.02 -0.05  0.11 -0.09
## PACF -0.03 -0.07 -0.04 -0.04 -0.09 -0.06  0.05 -0.03 -0.03  0.03  0.10 -0.01
acf2
## function (series, max.lag = NULL, plot = TRUE, main = NULL, ylim = NULL, 
##     na.action = na.pass, ...) 
## {
##     acf = stats::acf
##     pacf = stats::pacf
##     frequency = stats::frequency
##     num = length(series)
##     xfreq = frequency(series)
##     if (num > 59 & is.null(max.lag)) 
##         max.lag = max(ceiling(10 + sqrt(num)), 4 * xfreq)
##     if (num < 60 & is.null(max.lag)) 
##         max.lag = floor(5 * log10(num + 5))
##     if (max.lag > (num - 2)) 
##         stop("Number of lags exceeds number of observations")
##     if (is.null(main)) 
##         main = paste("Series: ", deparse(substitute(series)))
##     ACF = acf(series, max.lag, plot = FALSE, na.action = na.action, 
##         ...)$acf[-1]
##     PACF = pacf(series, max.lag, plot = FALSE, na.action = na.action, 
##         ...)$acf
##     LAG = (1:max.lag)/xfreq
##     if (plot) {
##         abline = graphics::abline
##         par = graphics::par
##         U = (-1/num) + (2/sqrt(num))
##         L = (-1/num) - (2/sqrt(num))
##         old.par <- par(no.readonly = TRUE)
##         if (is.null(ylim)) {
##             minA = min(ACF)
##             maxA = max(ACF)
##             minP = min(PACF)
##             maxP = max(PACF)
##             minu = min(minA, minP, L) - 0.01
##             maxu = min(max(maxA + 0.1, maxP + 0.1), 1)
##             ylim = c(minu, maxu)
##         }
##         par(mfrow = c(2, 1), cex.main = 1)
##         Xlab = ifelse(xfreq > 1, paste("LAG ÷", xfreq), "LAG")
##         tsplot(LAG, ACF, ylim = ylim, main = main, xlab = Xlab, 
##             ylab = "ACF", type = "h", ...)
##         abline(h = c(0, L, U), lty = c(1, 2, 2), col = c(8, 4, 
##             4))
##         tsplot(LAG, PACF, ylim = ylim, main = NULL, xlab = Xlab, 
##             ylab = "PACF", type = "h", ...)
##         abline(h = c(0, L, U), lty = c(1, 2, 2), col = c(8, 4, 
##             4))
##         on.exit(par(old.par))
##         ACF <- round(ACF, 2)
##         PACF <- round(PACF, 2)
##         return(rbind(ACF, PACF))
##     }
##     else {
##         return(cbind(ACF, PACF))
##     }
## }
## <bytecode: 0x7fbdc75b3d78>
## <environment: namespace:astsa>

En este primer modelo M1 se obtiene que no hay ruido blanco por lo que se deben explorar otras opciones para ajustar el modelo.

Se usa la función autoarima para obtener aproximación de parámetros en la modelación:

auto.arima(muestra)
## Series: muestra 
## ARIMA(2,0,1)(0,0,1)[12] with non-zero mean 
## 
## Coefficients:
##           ar1     ar2     ma1     sma1    mean
##       -0.5299  0.2388  0.7587  -0.8854  0.1154
## s.e.   0.1813  0.0498  0.1815   0.0291  0.0195
## 
## sigma^2 = 5.337:  log likelihood = -1191.2
## AIC=2394.4   AICc=2394.56   BIC=2419.98
M2<-arima(muestra,order=c(2,0,1),seasonal=list(order=c(0,1,1),period=12))
M2
## 
## Call:
## arima(x = muestra, order = c(2, 0, 1), seasonal = list(order = c(0, 1, 1), period = 12))
## 
## Coefficients:
##          ar1      ar2      ma1     sma1
##       0.6416  -0.0587  -0.4695  -1.0000
## s.e.  0.7543   0.1486   0.7525   0.0172
## 
## sigma^2 estimated as 9.132:  log likelihood = -1317.93,  aic = 2645.85
bic=AIC(M2,k = log(length(muestra)))
bic
## [1] 2667.17
acf2(M2$residuals)

##      [,1] [,2]  [,3] [,4]  [,5] [,6]  [,7]  [,8] [,9] [,10] [,11] [,12] [,13]
## ACF     0    0 -0.02 0.05 -0.03 0.03 -0.01 -0.05    0  0.07     0 -0.43  0.04
## PACF    0    0 -0.02 0.05 -0.03 0.03 -0.01 -0.05    0  0.07     0 -0.43  0.06
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF  -0.01  0.08  0.05  0.03 -0.03  0.03  0.06 -0.03 -0.07  0.04 -0.06 -0.02
## PACF  0.01  0.08  0.10 -0.01  0.00  0.01  0.03 -0.04 -0.02  0.06 -0.31  0.04
##      [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37]
## ACF  -0.01 -0.05 -0.09  0.00 -0.01 -0.01 -0.01  0.08  0.07 -0.09  0.08 -0.06
## PACF -0.02  0.00 -0.01 -0.01 -0.02 -0.01  0.06  0.04  0.08 -0.07 -0.13 -0.04
##      [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## ACF  -0.02  0.01  0.05 -0.02  0.07  0.05 -0.04 -0.05 -0.01  0.09 -0.03
## PACF -0.04  0.03  0.02  0.02  0.06  0.06 -0.01  0.02  0.06  0.03 -0.07
acf2
## function (series, max.lag = NULL, plot = TRUE, main = NULL, ylim = NULL, 
##     na.action = na.pass, ...) 
## {
##     acf = stats::acf
##     pacf = stats::pacf
##     frequency = stats::frequency
##     num = length(series)
##     xfreq = frequency(series)
##     if (num > 59 & is.null(max.lag)) 
##         max.lag = max(ceiling(10 + sqrt(num)), 4 * xfreq)
##     if (num < 60 & is.null(max.lag)) 
##         max.lag = floor(5 * log10(num + 5))
##     if (max.lag > (num - 2)) 
##         stop("Number of lags exceeds number of observations")
##     if (is.null(main)) 
##         main = paste("Series: ", deparse(substitute(series)))
##     ACF = acf(series, max.lag, plot = FALSE, na.action = na.action, 
##         ...)$acf[-1]
##     PACF = pacf(series, max.lag, plot = FALSE, na.action = na.action, 
##         ...)$acf
##     LAG = (1:max.lag)/xfreq
##     if (plot) {
##         abline = graphics::abline
##         par = graphics::par
##         U = (-1/num) + (2/sqrt(num))
##         L = (-1/num) - (2/sqrt(num))
##         old.par <- par(no.readonly = TRUE)
##         if (is.null(ylim)) {
##             minA = min(ACF)
##             maxA = max(ACF)
##             minP = min(PACF)
##             maxP = max(PACF)
##             minu = min(minA, minP, L) - 0.01
##             maxu = min(max(maxA + 0.1, maxP + 0.1), 1)
##             ylim = c(minu, maxu)
##         }
##         par(mfrow = c(2, 1), cex.main = 1)
##         Xlab = ifelse(xfreq > 1, paste("LAG ÷", xfreq), "LAG")
##         tsplot(LAG, ACF, ylim = ylim, main = main, xlab = Xlab, 
##             ylab = "ACF", type = "h", ...)
##         abline(h = c(0, L, U), lty = c(1, 2, 2), col = c(8, 4, 
##             4))
##         tsplot(LAG, PACF, ylim = ylim, main = NULL, xlab = Xlab, 
##             ylab = "PACF", type = "h", ...)
##         abline(h = c(0, L, U), lty = c(1, 2, 2), col = c(8, 4, 
##             4))
##         on.exit(par(old.par))
##         ACF <- round(ACF, 2)
##         PACF <- round(PACF, 2)
##         return(rbind(ACF, PACF))
##     }
##     else {
##         return(cbind(ACF, PACF))
##     }
## }
## <bytecode: 0x7fbdc75b3d78>
## <environment: namespace:astsa>
M3<-sarima(muestra,p=2,d=0,q=1,P=0,D=1,Q=4,S=12)
## initial  value 1.636228 
## iter   2 value 1.300917
## iter   3 value 1.138853
## iter   4 value 1.118309
## iter   5 value 1.102626
## iter   6 value 1.081712
## iter   7 value 1.074912
## iter   8 value 1.068452
## iter   9 value 1.065561
## iter  10 value 1.064848
## iter  11 value 1.063966
## iter  12 value 1.063702
## iter  13 value 1.063667
## iter  14 value 1.063600
## iter  15 value 1.063399
## iter  16 value 1.062820
## iter  17 value 1.062453
## iter  18 value 1.061801
## iter  19 value 1.061584
## iter  20 value 1.061167
## iter  21 value 1.061067
## iter  22 value 1.060661
## iter  23 value 1.060520
## iter  24 value 1.060175
## iter  25 value 1.060071
## iter  26 value 1.060058
## iter  27 value 1.060014
## iter  28 value 1.059982
## iter  29 value 1.059933
## iter  30 value 1.059915
## iter  31 value 1.059911
## iter  32 value 1.059911
## iter  32 value 1.059911
## iter  32 value 1.059911
## final  value 1.059911 
## converged
## initial  value 1.028001 
## iter   2 value 1.009554
## iter   3 value 1.006791
## iter   4 value 0.986489
## iter   5 value 0.981427
## iter   6 value 0.949323
## iter   7 value 0.946465
## iter   8 value 0.945083
## iter   9 value 0.940426
## iter  10 value 0.939957
## iter  11 value 0.939227
## iter  12 value 0.938873
## iter  13 value 0.938389
## iter  14 value 0.938212
## iter  15 value 0.938199
## iter  16 value 0.938198
## iter  17 value 0.938197
## iter  18 value 0.938194
## iter  18 value 0.938194
## final  value 0.938194 
## converged

M3$ttable
##          Estimate     SE  t.value p.value
## ar1        0.6895 1.1201   0.6156  0.5384
## ar2       -0.0455 0.2721  -0.1671  0.8674
## ma1       -0.4902 1.1208  -0.4374  0.6620
## sma1      -1.8247 0.1530 -11.9232  0.0000
## sma2       0.7507 0.1551   4.8415  0.0000
## sma3       0.0899 0.0998   0.9010  0.3680
## sma4      -0.0137 0.0488  -0.2817  0.7783
## constant   0.0004 0.0002   2.1393  0.0329
M3$AIC
## [1] 4.749353
M3$BIC
## [1] 4.823744

Para verificar si los residuos del modelo son ruido blanco, es decir validar si no hay patrones de autocorrelación se aplica la prueba Ljung-Box, con el fin de indentificar si la serie de tiempo exhibe autocorrelación significativa en varios lags. En este caso, los términos sma1 y sma2 son significativos, lo que sugiere que tienen un impacto estadísticamente significativo en el modelo.

Se tiene que el modelo M3 presenta buen ajuste, los coeficientes del modelo según la tabla t son significativos y los residuos del mismo se comportan como ruido blanco y presenta el mejor AIC.

CHEQUEO Y SELECCIÓN DEL MODELO:

Teniendo en cuenta que La hipótesis nula (H0) en la prueba de Box-Ljung es que no hay autocorrelación significativa en los residuos de la serie de tiempo hasta un cierto número de lags y La hipótesis alternativa (H1) es que hay autocorrelación significativa en los residuos hasta un cierto número de lags. La prueba calcula una estadística de prueba, conocida como estadístico Q de Ljung-Box, que se basa en los valores autocorrelacionados de los residuos hasta un número específico de lags. Se tiene que el modelo M4 presenta buen ajuste para la serie de tiempo ya que o hay autocorrelación significativa en los residuos de la serie de tiempo hasta un cierto número de lags y estos son ruido blanco.

Ahora bien alcomparar el BIC (Criterio de Información Bayesiano) y el AIC (Criterio de Información de Akaike) favorecen modelos más pequeños y parsimoniosos pero el BIC tiende a penalizar la complejidad más fuertemente que el AIC de los dos modelos se selecciona el modelo M3 por tener mejor BIC ``

Se puede apreciar gráficamente que el modelo M3 presenta un comportamiento similar a los datos de temperatura del aire, es decir que representa la serie de tiempo de temperatura media del aire de los últimos 44 años a una atmósfera de presión en el meridiando Greenwich, latitud cero.

Se aprecia estacionalidad año tras año si mayores fluctuaciones en los rendimientos año tras año teniendo en cuenta que los comportamientos mensuales a través del tiempo son bastante estables por lo tanto no se contemplará modelar la volatilidad

M4<-sarima(muestra, p=2,d=0,q=1,P=0,D=1,Q=4,S=12)
## initial  value 1.636228 
## iter   2 value 1.300917
## iter   3 value 1.138853
## iter   4 value 1.118309
## iter   5 value 1.102626
## iter   6 value 1.081712
## iter   7 value 1.074912
## iter   8 value 1.068452
## iter   9 value 1.065561
## iter  10 value 1.064848
## iter  11 value 1.063966
## iter  12 value 1.063702
## iter  13 value 1.063667
## iter  14 value 1.063600
## iter  15 value 1.063399
## iter  16 value 1.062820
## iter  17 value 1.062453
## iter  18 value 1.061801
## iter  19 value 1.061584
## iter  20 value 1.061167
## iter  21 value 1.061067
## iter  22 value 1.060661
## iter  23 value 1.060520
## iter  24 value 1.060175
## iter  25 value 1.060071
## iter  26 value 1.060058
## iter  27 value 1.060014
## iter  28 value 1.059982
## iter  29 value 1.059933
## iter  30 value 1.059915
## iter  31 value 1.059911
## iter  32 value 1.059911
## iter  32 value 1.059911
## iter  32 value 1.059911
## final  value 1.059911 
## converged
## initial  value 1.028001 
## iter   2 value 1.009554
## iter   3 value 1.006791
## iter   4 value 0.986489
## iter   5 value 0.981427
## iter   6 value 0.949323
## iter   7 value 0.946465
## iter   8 value 0.945083
## iter   9 value 0.940426
## iter  10 value 0.939957
## iter  11 value 0.939227
## iter  12 value 0.938873
## iter  13 value 0.938389
## iter  14 value 0.938212
## iter  15 value 0.938199
## iter  16 value 0.938198
## iter  17 value 0.938197
## iter  18 value 0.938194
## iter  18 value 0.938194
## final  value 0.938194 
## converged

M4$ttable
##          Estimate     SE  t.value p.value
## ar1        0.6895 1.1201   0.6156  0.5384
## ar2       -0.0455 0.2721  -0.1671  0.8674
## ma1       -0.4902 1.1208  -0.4374  0.6620
## sma1      -1.8247 0.1530 -11.9232  0.0000
## sma2       0.7507 0.1551   4.8415  0.0000
## sma3       0.0899 0.0998   0.9010  0.3680
## sma4      -0.0137 0.0488  -0.2817  0.7783
## constant   0.0004 0.0002   2.1393  0.0329
residuos4<-M4$fit$residuals
acf2(residuos4)

##      [,1] [,2]  [,3] [,4]  [,5] [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13]
## ACF     0 0.01 -0.03 0.04 -0.03 0.04 -0.01 -0.02 -0.02  0.09  0.01 -0.01  0.03
## PACF    0 0.01 -0.03 0.04 -0.03 0.04  0.00 -0.02 -0.01  0.08  0.01 -0.01  0.04
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF  -0.01  0.04  0.03     0 -0.03  0.02  0.04 -0.02 -0.01  0.03  0.01     0
## PACF -0.01  0.05  0.03     0 -0.02  0.02  0.04 -0.03  0.00  0.02  0.01     0
##      [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37]
## ACF  -0.01 -0.05 -0.07  0.00 -0.01     0 -0.02  0.06  0.07 -0.06  0.06 -0.01
## PACF -0.02 -0.05 -0.06 -0.01 -0.02     0 -0.02  0.05  0.08 -0.07  0.07  0.00
##      [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## ACF  -0.02 -0.02  0.03 -0.01  0.02  0.02 -0.06 -0.03  0.00  0.06 -0.02
## PACF -0.02 -0.02  0.02  0.00  0.03  0.02 -0.08 -0.01 -0.01  0.05 -0.01
Box.test(residuos4, lag=8, type='Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  residuos4
## X-squared = 3.0098, df = 8, p-value = 0.9337
Box.test(residuos4, lag=16, type='Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  residuos4
## X-squared = 9.5372, df = 16, p-value = 0.8897
Box.test(residuos4, lag=24, type='Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  residuos4
## X-squared = 11.871, df = 24, p-value = 0.9813

según estos resultados del test de Ljung-Box, no hay evidencia significativa de autocorrelación en los residuos del modelo SARIMA hasta los rezagos considerados (8, 16, y 24). Esto sugiere que el modelo SARIMA ha capturado adecuadamente la estructura de autocorrelación en la serie de tiempo. Lo anterior implica que se ha encontrado el Modelo.

Pronóstico

A continuación se presenta el resultado del análisis de la serie de tiempo de temperatura del aire, se identificó que la serie tiene un comportamiento estacional inicialmente no estacionaria que tuvo que diferenciarse estacionalmente para poder ser modelada, no fue necesario modelar volatilidad ni cambios estructurales en la serie y con la modelación se pudo hallar un modelo SARIMA con un componente autoregresivo (AR), un componente de media móvil (MA), tanto en la parte no estacional como en la estacional, y también incluye diferenciación para manejar tendencias no estacionarias y estacionalidad mensual. La elección específica de estos parámetros se basa en el análisis y la comprensión de la serie de tiempo particular que se ha modelado.

A continuación se presenta la gráfica del modelo con los datos de muestra frente a la serie original de temperatura del aire.

sarima.for(muestra,n.ahead=24,p=2,d=0,q=1,P=0,D=1,Q=4,S=12)
## $pred
##          Jan       Feb       Mar       Apr       May       Jun       Jul
## 45                                                                      
## 46 0.7767042 4.1919155 3.1136137 5.9711904 3.6871036 2.0276635 1.8500717
## 47 0.3802075 0.7470753 0.6233433 0.7555095 0.5028296 0.3187339 0.2930456
##          Aug       Sep       Oct       Nov       Dec
## 45                     1.4909667 3.6729078 2.0826290
## 46 1.2750707 0.5939942 0.6005659 0.7797546 0.5516416
## 47 0.2752712 0.2099281                              
## 
## $se
##         Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 45                                                                        
## 46 2.387892 2.389150 2.389627 2.389807 2.389876 2.389901 2.389911 2.389914
## 47 3.065040 3.065670 3.065909 3.066000 3.066034 3.066047 3.066052 3.066054
##         Sep      Oct      Nov      Dec
## 45          2.330450 2.376299 2.385897
## 46 2.389914 3.035797 3.058569 3.063387
## 47 3.066055
lines(Vwind.ts)

fig <- plot_ly(y = Vwind.ts, name = "Original Data", type = 'scatter', mode = 'lines') 
fig <- fig %>% add_trace(y = Vwind.ts+M4$fit$residuals, name = "Predicted Values", connectgaps = TRUE)
fig

Se hace la comparación final de los datos originales vs los datos generados por el modelo.